REMARKS 

The present invention is a method for modeling fluid flows in a fractured 
multiplayer porous medium to simulate interactions between pressure and flow rate 
variations in a well through the medium. The method discretizes the fractured 
medium by a mesh pattern with fractured meshes centered on nodes at fracture 
intersections with each node being associated with a matrix volume. Flows are 
determined between each fracture mesh and the associated matrix volume in a 
pseudosteady state by using an image processing algorithm. The use of the image 
processing algorithm as part of the determination of flows is described in paragraph 
[0029] of the Substitute Specification. 

Submitted herewith is a new Declaration of Bernard Bourbiaux pursuant to 
37 C.F. R. §1.132. Mr. Bourbiaux's Declaration is substantively identical to his 
previous Declaration except he has substituted reference to French Patent 
2,757,947 which is prior art as noted by the Examiner in Section 15-1 . However, in 
Mr. Bourbiaux's Declaration it is noted that there is a typographical error in the 
quotation from the specification of the present application in that, under the heading 
"(a) Block geometry", there is a reference to French Patent 2,757,957 which should 
be French Patent 2,757,947 and furthermore, in the second paragraph, under the 
quotation from paragraph [0027], reference to United States Patent 6,023,056 
appears. That reference should have been replaced with the aforementioned 
French Patent 2,757,947. It is requested that the Examiner consider the references 
to United States Patent 6,023,056 and the typographical error regarding French 
Patent 2,757,957 to be references to French Patent 2,757,947 throughout in 
evaluating the merits of the Declaration. Applicants will submit a third Declaration 
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correcting these errors if the Examiner requires another declaration that this is 
necessary. 

Section 4 of the Office Action has objected to the Substitute Specification filed 
March 3, 2004 regarding dots in equations not being properly positioned. The 
Second Substitute Specification submitted herewith deletes dots in the equations. 
The references to United States patents, which the Examiner correctly notes are not 
prior art, have been replaced in the Second Substitute Specification with the 
counterpart French patents which are prior art and therefore, are information 
relevant to the issue of enablement raised by the Examiner that is discussed below. 

Claims 7-12 were rejected under 35 U.S.C. §112, first paragraph, as allegedly 
not being enabling. 

Submitted herewith is the Declaration of Bernard Bourbiaux pursuant to 
37 C.F.R. §1.132 which addresses the Examiner's conclusions in Section 7-1 of the 
Office Action that the methodology of pore volume calculation of fracture matrices 
has not been disclosed in the specification. Mr. Bourbiaux's Declaration is 
substantively identical to his previous Declaration except that reference to United 
States Patent 6,023,652 has been replaced with corresponding French Patent 
2,757,947 as discussed above. 

Mr. Bourbiaux's Declaration emphasizes relevant portions of the text of 
paragraphs [0020]-[0027] also appearing in the Second Substitute Specification 
submitted herewith which should be specifically considered by the Examiner in 
addition to the overall text of paragraphs [0020]-[0027] of the specification. 
Mr. Bourbiaux's Declaration concludes the calculation of volume as discussed in 
paragraphs [0020] and [0027] is a simple volumetric calculation of area times height 
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as understood by persons of ordinary skill in the field of the present invention at the 
time the application was filed. 

Mr. Bourbiaux's Declaration describes French Patent 2,757,947, which is 
referred to in the Second Substitute Specification of the present application, 
regarding its relevant teachings pertaining to the calculation of volume. 
Mr. Bourbiaux provides a diagram of fractured meshes, which is consistent with 
Fig. 2 of French Patent 2,757,957 on page 5 of his Declaration, which explains 
relevant dimensions which are used to compute the area of a fracture mesh volume. 
Mr. Bourbiaux concludes on page 5 that "[i]f h is the thickness of the geological 
layer containing the considered fracture mesh, and h and l 2 are the respective 
lengths (heights) of the two constitutive fracture segments A1B1 and A2B2, and ai 
and a 2 their respective apertures, then the fracture mesh volume is very close to..." 
which is the equation reproduced at the bottom of page 5 of his Declaration. 

It is submitted that based upon Mr. Bourbiaux having an engineering degree 
and a Masters degree and further having worked in research relating to modelling of 
fractured reservoirs and simulation thereof, that his Declaration is entitled to status 
as an expert opinion regarding the sufficiency of disclosure regarding calculation of 
volume. Moreover, as described above, it is submitted that Mr. Bourbiaux's 
Declaration demonstrates that a person of ordinary skill in the art, without undue 
experimentation, would be enabled to calculate volume without undue 
experimentation. 

Claims 9 and 12 stand rejected under 35 U.S.C. §112, second paragraph, as 
being indefinite for failing to particularly point out and distinctly claim the subject 
matter which Applicant regards as the invention. Claim 9 has been amended to 
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depend from claim 7 and recites each fractured layer as discretized in pixels and the 
matrix volume associated with each fracture mesh is defined in determining a 
distance from each pixel to a closest fracture mesh. It is submitted that this subject 
matter is definite and further limits claim 7. 

Claims 10 and 1 1 stand rejected under 35 U.S.C. §112, second paragraph, 
as being incomplete for omitting essential steps with the Examiner asserting that the 
essential step is "determining a distance". This ground of rejection is traversed for 
the following reasons. As the Examiner is aware, the second paragraph of 35 
U.S.C. §112 provides that the Applicant's claims should particularly point out and 
distinctly claim the subject matter which the applicant regards as the invention . It is 
submitted that the methodology of determining the distance, as asserted by the 
Examiner to be an essential step, is not necessary since there is no question that 
the specification supports how the determination of the distance is made. Applicant's 
position is that it is not necessary in dependent claims 10 and 1 1 to recite the 
methodology of determining the distance as asserted by the Examiner to be an 
essential step and would be unduly limiting. 

The Examiner should note that newly submitted claims 13 and 14 further limit 
claims 10 and 1 1 in reciting "discretizing each fracture layer and pixel; and 
determining a distance from each pixel to a closest fracture mesh and wherein the 
determining at any point is from the distance." These limitations recite the 
"determining a distance" which the Examiner asserts is a necessary step in 
claims 10 and 11. 

Claims 7 and 8 stand rejected under 35 U.S.C. §103 as being unpatentable 
over United States Patent 6,023,656 (Cacas et al) in view of United States Patent 
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6,094,619 (Noetinger et al). This ground of rejection is traversed for the following 
reasons. 

Claim 7 recites a method for modelling fluid flows in a fractured multilayered 
porous medium to simulated interactions between pressure and flow rate variations 
in a well through the medium, comprising: discretizing the fractured medium by a 
mesh pattern with fracture meshes centered on nodes at fracture intersections with 
each node being associated with a matrix volume; and determining flows between 
each fracture mesh and the associated matrix volume in a pseudo-steady state by 
using an image processing algorithm. Neither Cacas et al or Noetinger et al disclose 
the aforementioned "determining flows between each fracture mesh and the 
associated matrix volume in a pseudo-steady state by using an image processing 
algorithm." Noetinger uses a different mechanism which is the so-called random 
walk methodology which is described in column 6, lines 25-35, for determining the 
flows. In practice, the random walk methodology has turned out to be both a time 
consuming method, is hard to use with the discretized medium of the present 
invention and is distinct from the claimed image processing algorithm. Cacas et al 
do not disclose the utilization of the image processing algorithm for the claimed 
determination of flows. See the description of step 104 in Fig. 10A. Accordingly, if 
the proposed combination of Cacas et al and Noetinger et al were made, the subject 
matter of claim 7 would not be achieved. 

Claims 9-12 stand rejected as being unpatentable over Cacas et al, 
Noetinger et al and United States Patent 6,064,944 (Sarda et al). Sarda has been 
cited as disclosing a method based on a pixel representative which determines the 
dimensions of equivalent blocks. Sarda does not cure the deficiencies of 



10 



Cacas et al and Noetinger as discussed above. Moreover, it is submitted that the 
Examiner's proposed combination is based upon impermissible hindsight. 

In view of the foregoing amendments and remarks, it is submitted that each of 
the claims in the application is in condition for allowance. Accordingly, early 
allowance thereof is respectfully requested. 

To the extent necessary, Applicants petition for an extension of time under 
37 C.F.R. §1.136. Please charge any shortage in fees due in connection with the 
filing of this paper, including extension of time fees, to Deposit Account No. 01-2135 
(612.37806X00) and please credit any excess fees to such Deposit Account. 



Respectfully submitted, 




ANTONELLLTERRY, STOUT & KRAUS, LLP 



Donald E. Stout 
Registration No. 26,422 
(703)312-6600 
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ABSTRACT 



M e thod A method is disclosed for modelling fluid flows in a fractured multilayer 
porous medium by taking account o f accounting for the real geometry of the fracture 
network and the local exchanges between the porous matrix and the fractures at each 
node of the network, thus allowing simulation of the interactions between the pressure 
and flow rate variations in a well running across saidthe medium. The method 
essentially comprises discr e tization o f discretizing the fractured medium by means of a 
mesh pattern, with fracture meshes c e ntr e d centered on nodes at the various intersections 
of the fractures^ with each node being associated with a matrix volume, and 
determination of the flows between each fracture mesh and the associated matrix 
volume in a pseudosteady state. The method can be applied in hydrocarbon production 
well testing for e xampl e. 
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BACKGROUND OF THE INVENTION 
FIELD OF THE INVENTION Field of the Invention 

The present invention relates to a method for simulating well tests on a discrete 
fractured reservoir model. 

5 The method can b e implemented for exampl e in th e field of oil production by 

reservoir e ngin ee rs in ord e r to obtain r e liabl e flow prediction s . 

A w e ll test is to b e s imulat e d in a p e rm e able porous m e dium (r e s e rvoir) cross e d 
through by a n e twork of fractur es much mor e conducting than th e porous matrix. 
Fractur e d r e s e rvoirs constitut e an e xtr e m e type of h e t e rog e neous r e servoirs comprising 
10 two v e ry diff e rent m e dia, a matrix medium containing the most part of th e oil in plac e 
and having a low p e rmeab i lity, and a fractured m e dium r e presenting l e ss than 1 % of 
the oil in plac e and highly conducting. The fractured medium itsel f c an be complex, 
with diff e r e nt series of fractur e s charact e riz e d by th e ir r e sp e ctive d e nsity, l e ngth, 
ori e ntation, inclination and opening. 

15 During w e ll t e sting, the flow rat e conditions impos e d on the well lead th e oil 

contain e d in th e r e s e rvoir to flow toward s th e w e ll. It i s a singl e phase (the only mobil e 
phas e is th e oil) and compressibl e flow (th e rock of the res e rvoir and the oil fluid ar e 
compr e ssibl e ). 

For any e l e m e ntary volume of the r e s e rvoir, th e pressure of the oil contained in this 
20 vo l um e i s gov e rn e d by the following equation, if gravity is not taken into account : 



wh e r e- : 

§ : r e pr e sent s th e por e volum e , 

C+, the total compr es sibility (fluid + rock), 

K, th e p e rm e ability of the rock, 

u, the viscosity of th e fluid, 

Q, th e incoming flow, which is zero everywh e r e e xcept in plac e s where th e we l l 
communicate s with th e reservoir, and 
P, the unknown pressur e . 

In ord e r to simulat e a w e ll t e st, what e v e r the medium , this e quation ha s to b e 
s olved in spac e and in tim e . Discretization of the reservoir (m e sh patt e rn) is therefor e 
p e rform e d and solution of th e problem consists in finding the pressur es of th e m e sh es 
with tim e , itself di s cr e tized in a certain number of tim e int e rvals. 

In the ca s e of a fractured medium, if all the compl e xity of the fracture network is 
to be tak e n into account, it is nec e s s ary to have a mesh patt e rn that r e presents this 
network accurately. Furth e rmore, in such a medium, th e p e rm e abiliti e s K are g e n e rally 
much high e r in the fractures than in th e matrix, and th e flow s are cons e quently fa s t in 
th e fractur e s and s low in th e matrix. 

BACKGROUND OF THE INVENTlON Description of the Prior Art 

Single-phase flow computing tools are currently used ; howeve r . However , 
the vthese tools are not applied to the real geologic medium in aM-4t sview of the 
complexity bt rtthereof and instead are applied to a homogenized representation, 
according to the reservoir model called a "double-medium model" described for 
example by Warren and Root in "The Behavior of Naturally Fractured Reservoirs", SPE 




Journal, September 1963. Any elementary volume of the fractured reservoir is thus 
modelled in the form of a set of identical parallelepipedic blocks limited by an 
orthogonal system of continuous uniform fractures oriented in the direction of one of 
the three main directions of flow (model referred to as "sugar box" model). The fluid 
5 flow on the reservoir scale occurs through the fractured medium only, and fluid 
exchanges take place locally between the fractures,, and th e The matrix blocks are not 
applied to the complex medium. This representation, which does not reproduce the 
complexity of the fracture network in a reservoir, is however effective but at the level of 
a reservoir mesh whose typical dimensions are 100 m x 100 m. 

10 It is however preferable to carry out the flow calculations on the "real" geologic 

model instead of an equivalent homogeneous model, which affords the double 
advantage of allowing : 

to validate validation of the geologic image of the reservoir bwt tprovided by the 

geologist from all the information he has gathered on the reservoir fracturation (this 
15 validation being performed by comparison with the well test data), 

te reliably testtesting the sensitivity of the hydraulic b e haviou r behavior of the 

medium to the uncertainties on the geologic image of the fractured medium. 

A well-known modelling method con s i s ts in| s finely meshing the fracture network 
and the matrix while making no approximation concerning fluid exchanges between the 
20 two media. It is however difficult to implement because the often complex geometry of 
the spaces between the fractures makes it difficult to me s h th e m and. the meshing of the 
spaces, mln any case, the number of meshes to be processed is finally often 
trem e ndous very large . The complexity increases still further with a 3D mesh pattern. 
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Techniques for modelling fractured porous media are described in US patents ¥Rr- 
A 2,757.947 6.023.656 and FR A 2.757.957 6.064.944 filed by the apptieant assignee . 
The form e r technique of the first patent relates to determination of the equivalent 
fracture permeability of a fracture network in an underground multilayer medium from 

5 a known representation of thisthe network, allowing te — sy s tematically 
connect svstematic connection of fractured reservoir characterization models to double- 
porosity simulators in order to obtain more realistic modelling of a fractured 
underground geologic structure. The s e cond technique of the second patent relates to 
simplified modelling of a porous heterogeneous geologic medium (such as a reservoir 

10 crossed through by an irregular network of fractures for example) in the form of a 
transposed or equivalent medium so that the transposed medium is equivalent to the 
original medium, according to a determined type of physical transfer function (known 
for the transposed medium). 

SUMMARY OF THE INVENTION 

15 The method according to the invention allows to simulat e simulation single-phase 

flows in a fractured medium. 

The method can be implemented for example in the field of oil production by 
reservoir engineers in order to obtain reliable flow predictions. 

A well test is to be simulated in a permeable porous medium (reservoir) crossed 
20 through by a network of fractures much more conducting than the porous matrix. 
Fractured reservoirs constitute an extreme type of heterogeneous reservoirs comprising 
two very different media, a matrix medium containing the most part of the oil in place 
and having a low permeability, and a fractured medium representing less than 1 % of 
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the oil in place and highly conducting. The fractured medium itself can be complex, 
with different series of fractures characterized bv their respective density, length, 
orientation, inclination and opening. 

During well testing, the flow rate conditions imposed on the well lead the oil 
5 contained in the reservoir to flow towards the well. It is a single-phase (the only mobile 
phase is the oil) and compressible flow (the rock of the reservoir and the oil fluid are 
compressible). 

For any elementary volume of the reservoir, the pressure of the oil contained in this 
volume is governed bv the following equation, if gravity is not taken into account : 

bPjS^- = div(—J7P) + Q (1) 
o t p. 

10 where : 

(b : represents the pore volume, 
C T . the total compressibility (fluid + rock), 
K. the permeability of the rock, 
u, the viscosity of the fluid. 
15 O, the incoming flow, which is zero everywhere except in places where the well 
communicates with the reservoir, and 
P. the unknown pressure. 

In order to simulate a well test, whatever the medium, this equation has to be solved 
in space and in time. Defining of the reservoir (mesh pattern) is therefore performed and 
20 solution of the problem finds the pressures of the meshes with time, itself defined in a 
certain number of time intervals. 




In the case of a fractured medium, if all the complexity of the fracture network is to 
be taken into account it is necessary to have a mesh pattern that represents this network 
accurately. Furthermore, in such a medium, the permeabilities K are generally much 
higher in the fractures than in the matrix, and the flows are consequently fast in the 
5 fractures and slow in the matrix. 

The obj e ct of th e method according to the invention i s to mod e l models fluid flows 
in a fractured multilayer porous medium by taking account o f accounting for the real 
geometry of the fracture network and the local exchanges between the porous matrix 
and the fractures at each node of the network, and thus to allow simulation o f simulates 

10 the interactions between the pressure and flow rate variations in a well running across 
sa4dthe medium. It is characterized in that th e The fractured medium is discretized by 
means of a mesh pattern wherein the fracture meshes are c e ntr e d centered on nodes at 
the various intersections of the fractures, each node being associated with a matrix 
volume, and the flows between each fracture mesh and the associated matrix volume are 

15 determined in a pseudosteady state. 

The matrix volume associated with each fracture mesh in each layer is for example 
delimited defined by all of the points that are closer to the corresponding node than to 
n e ighbourin gn eighboring nodes. 

In order to determine the matrix volume associated with each fracture mesh, each 
20 fractured layer is for example discretized in pixels and the distance from each pixel to 
the closest fracture mesh is determined. 
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The value of the transmissivity for each fracture mesh - matrix block pair is 
determined for example by considering that the pressure varies linearly as a function of 
the distance from the point considered to the fracture mesh associated with the block. 

The method according to the invention greatly simplifies calculation of the 
5 exchanges between the matrix and the fractures. In fact, there is no specific mesh 
pattern for the matrix, but each node of the fracture network is associated with a matrix 
volume, and the flow between each fracture node and its associated block is calculated 
in a pseudosteady state. The innovation invention concerns calculation of the volume of 
each block and of thea proportionality coefficient (called transmissivity) that connects 
10 the matrix-fracture flow rate to the matrix-fracture pressure difference. 

The main advantage of the method according to the invention in relation to methods 
applied to a finely meshed real fracture network lies, as detailed in the description 
hereafter, in thea considerable reduction in thea number of meshes c e ntr e d centered on 
fracture nodes used to solve equation (1). It can be checked that th e The simulation rate 
15 saving in relation to prior fine-mesh methods i scan be confirmed to be greater than the 
simple arithmetical ratio of the number of meshes used in both cases. 

BRIEF DESCRIPTION OF THE DRAWINGS 

Other features and advantages of the method according to the invention will be 
clear from reading the description hereafter of a non limitative realisation example, with 
20 reference to the accompanying drawings wherein : 

- EigHf eFig. 1 shows an example of a 2D fracture mesh, FM, 

- EigufeFig. 2 shows an example of a 2D matrix block, MB, 




£igufe Fig, 3 shows an imposed flow rate variation curve F(t) in a well test, and 

Ftgttfe Fig, 4 shows an example of a fracture network for which the method 
according to the invention leads to a radiea lsubstantiat reduction in the number of 
meshes to be processe d, and 

5 - Figs. 5 and 6 illustrate flow charts of the method of the present invention . 

DETAILED DESCRIPTION OF THE M ETHOD INVENTION 
1) Fracture network discretization 

The method comprises meshing the multilayer fractured reservoir modelled for 
example by means of the fractured reservoir modelling technique described notably in 
10 the aforementioned US patent FR 2,757,947 6.023.656 , assuming that the fractures are 
substantially perpendicular to the layer boundaries. 

In order to form this mesh pattern, all the fracture intersections or computing nodes 
GN are located (Fig.l) and fracture meshes are formed by associating a single matrix 
block MB (Fig.2) with each node. For the purpose of 3D modelling, additional nodes 

1 5 have to be added in the neighbourin g neighboring layers in order to take account effor 
the fractures that run across several of th e m layers. The c om p utin g computed nodes thus 
defined constitute the eentf ecenter of the fracture meshes £M. The meshes are given 
boundaries such as the ends of the fractures on the one hand and the midpoint of the 
segments connecting the computing nodes on the other hand. Considering these 

20 boundaries imposed on the meshes, thewthe volume <j> of the meshes can be calculated 
(seeSee relation 1). Calculation of the connections between fracture meshes 
(transmissivities) used for calculation of the inter-mesh flows can be carried out 
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according to the method described in the aforementioned US patent FR- 
2,757.91 7 6,023.656 . 

2) Matrix medium discretization 

Discretization of the matrix medium consist s in assignin g assigns a rock volume to 
5 each fracture mesh. During dynamic simulation, exchanges between the matrix and the 
fractures wU4-beare calculated in the pseudosteady state (exchange flow proportional to 
the pressure difference) between each fracture mesh and itsthe associated single matrix 
block. 

For calculation of thgsethe matrix volumes, the problem is dealt with layer by layer, 
10 which amounts to disregarding the vertical flows in relation to the horizontal flows in 
the matrix. 

For a given layer, the blocks are defined as follows : 

the block heights are equal and fixed to the height of the layer, 

in the plane of the layer, the surface area of the block associated with a fracture 
1 5 mesh is all the points that are closer to thisthe fracture mesh than to another one. 

Physically, this definition implies that the boundaries at zero flow in the matrix are 
the equidistance lines between the fracture meshes. This approximation is valid if the 
n e ighbourin g neighboring fractures are at very close pressures, which is true in the 
presence of highly conducting fractures in relation to the matrix. 

20 a) Block geometry 

In order to determine the geometry of the blocks in a given layer, a two-dimensional 
problem consisting in finding has to solved which finds , for each fracture mesh, the 
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points that are closer to this mesh than to another one has to b e solv e d . This problem is 
solved by using for example, but preferably, the geometric method described in the 
other aforementioned US patent FR 2,757.957 6,064,944 . which allows, by discretizing 
the fractured layer into a set of pixels and by applying an image processing algorithm, 
5 to determine the distance from each pixel to the closest fracture. During the 
initialization phase of this algorithm, value 0 (zero distance) is assigned to the pixels 
belonging to a fracture and a high value is assigned to the others. Furthermore, if the 
number of the fracture mesh to which each fracture pixel belongs is given in this phase, 
the same algorithm finally allows to determine, for each pixel : 

10 - the distance from this pixel to the closest fracture mesh, 

the number of the closest fracture mesh. 

All of the pixels having the same fracture mesh number constitute the surface area 
of the matrix block associated with this mesh. Multiplying this surface area by the 
height of the layer allows to obtain obtaining the volume of the block associated with the 
15 mesh. Figure Fig. 2 shows an example of a thus obtained 2D matrix block. 

By construction, the sum of the surface areas of the blocks is equal to the total 
surface area of the layer, which guarantees that no volume is added to or taken from the 
layer. 

For each matrix block, the image processing algorithm also gives the distance from 
20 each pixel of the block to the associated fracture mesh. This data is used to calculate the 
transmissivity between the fracture mesh and the matrix block. 
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b) Matrix-fracture transmissivity 

On the assumption of a pseudosteady state where the exchange flow is considered 
to be proportional to the difference in the pressures of each one of them, we have the 
following relation exists : 



F mf is the matrix-fracture flow, 

T mf , the matrix-fracture transmissivity, 

jx, the viscosity, 

P m , the pressure of the matrix block, and 
10 P f , the pressure of the associated fracture mesh. 

On the assumption of a pseudosteady state, the value of transmissivity T mf is 
constant during simulation. This value is here calculated for each fracture mesh - matrix 
block pair. 

For this calculation, it is assumed that, in the matrix block, the pressure varies 
15 linearly as a function of the distance from the point considered to the fracture mesh 
associated with the block. The transmissivity is defined as follows : 
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where : 
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D 



where : 




I is the length of the fractures in the fracture mesh, 
H, the height of the layer (and of the matrix block), 
K, the permeability of the matrix, and 

D, the distance from the fractures for which the pressure is the mean pressure of the 
5 block. 

1, H and K are known (factor 2 comes from the fact that the fracture has two faces). 
Calculation of D is made from the information provided by the image processing 
algorithm concerning the distances from the pixels to the closest fractures. In fact, 
according to the hypothesis of linear variation of the pressure as a function of the 
10 distance to the fracture, w e hav e the relationship is : 

D - -Jd(s\ds (4) 

s s 

where S is, in 2D, the surface area of the matrix block and d(s) is the function giving the 
distances between the points of the matrix block and the associated fracture mesh. 

In terms of pixels, w e thus s imply have the relationship is : 

D = jr^ dn 44 ^— 

15 where N is the number of pixels of the matrix block and d n the distance from pixel n to 
the fracture mesh. 

3) Well representation 

A well is represented by its geometry on the one hand and by the flow rate 
constraints imposed thereon with time on the other hand. 
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a) Geometry 

A well consists ofl s a series of connected segments that intersect the network 
fractures. The geometric representation of a well is therefore a 3D broken line. Some 
segments of this broken line communicate with the reservoir. The points of 
5 communication between the well and the network are the points of intersection of these 
segments communicating with the fractures. 

b) Flow rate constraints 

The flow rate constraints are imposed at the point where the well enters the 
fractured medium. They The flow rates are entered by the user in the form of a step 
10 function giving the flow rate as a function of time. The flow rate change times of the 
well indicate the various periods to be simulated. 

For a conventional well test (Fig. 3), a flow rate F(t) is imposed for a first period 
after which the well is closed (zero flow rate) and the pressure buildup is observed in 
the well. 

15 4) Dynamic simulation 

During dynamic simulation, the simulated period of time is split up into time 
intervals whose duration ranges, for a well test, between 1 s (just after opening or 
closing of the well) and sem ein a range of hours. 

Passage from the time n (where all the pressure values are known) to the time n+ 1 
20 is achieved by solving, in all the meshes and blocks of the reservoir, the following 
equation which corresponds to relation 1 once discretized : 




_ pn p 

where : 

i, the number of the mesh or of the block, 
j, the number of the mesh or of the block next to i, 
5 P n , the pressure at the time n, 
t n , the time at the time n, 
Qi, the flow entering i, and 
Ty, the connection transmissivity between i and j. 

Apart from the pressures, the terms of this equation are known. The pore volumes 
10 of the fracture meshes and of the matrix blocks (<|>j) are known by means of the mesh 
pattern, the fracture-fracture and matrix-fracture transmissivities (Ty) are calculated as 
described above and flow rates (Q s ) are zero everywhere except at the well-reservoir 
connections where they are imposed. The mesh pattern is known using the method 
described in the above cited US patent 6.023.656. Fracture lengths and openings 
15 (widths) are given by the characteristics of the fracture network and are input data. The 
pore volume in a given fracture mesh is determined from the length of fractures within 
the mesh and the corresponding opening of the fracture. 

a) Elimination of the unknowns relative to the matrix 

Equation (6)(7) can be solved by inversion of a linear system where the unknowns 
20 are the pressures of the fracture meshes and the pressures of the matrix blocks at the 
time n+ 1 . 
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However, since each matrix mesh is only connected to one matrix block and since 
the blocks are not connected to each other, it is possible to eliminate the unknowns of 
the blocks of the system whose size is then divided by 2. 

In fact, for a matrix block, equation (Q(7) is expressed as follows : 

I [ 

5 subscripts m and f denoting the matrix and the fracture respectively, and C Tm being the 
total compressibility of the matrix. 

We hav e The relationship is Cthi — c m + Sw.c w + (1-Sw).c 0 , where Sw is the residual 
water saturation in the matrix, c m the compressibility of the matrix, c w the 
compressibility of the water and c 0 the compressibility of the oil, and § m is the pore 
1 0 volume of the matrix block, irer that is the volume of the block multiplied by the porosity 
of the matrix. 

We The relationship is deduce therefrom that : 

t 1 = , \ T j >: + , Um rr .j > n + i (8) 



with : 



T 

Itt 

By transferring this expression into the equation relative to the fracture mesh, we 
15 ebtai ftthe relationship is obtained : 
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J J j r 1 



m 



m 



m 



(9) 



where i is the number of the fracture mesh considered, j the numbers of the fracture 
meshes connected to i, and m the number of the block associated with fracture mesh i. 

Furthermore, we have the relationship is : 



5 where C T f is the total fracture compressibility : C Tf = c f + c 0 . 

Reduction of the number of unknowns by a factor 2 is a further advantage of the 
present method. The solution time for a linear system develops approximately like the 
square of the number of unknowns, and consequently elimination of the matrix 
unknowns leads to a considerable saving in the rate of simulation of a well test on a 
10 fractured network. 

b) Time interval solution algorithm 

The general algorithm for solving a time interval con s ist s in formin gforms a linear 
system corresponding to equation (9) relative to any fracture mesh i, in solving this 
linear system and in updating the unknowns. 

1 5 The stages of this algorithm are as follows : 

Forming the linear system with calculation of the terms of equation (9) 
Calculation of the fracture medium contribution 
accumulation term (A ; ) 




n 
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connection between fracture meshes (Ty/u) 



well boundary conditions (Qi) 



Calculation of the matrix medium contribution 



accumulation term (A m ) 



5 



matrix-fracture connection term (U m ) 



Solving the linear system (so as to find the fracture pressure values at the time n+1) 
by means of the iterative conjugate gradient algorithm (CGS) 

- Updating with calculation of the matrix block pressures using equation (8), and time 
interval management. 



During well test simulation, the time intervals are generally very short after opening 
or closing of the well, because the pressure variations are high at these times. Later, the 
time intervals can be longer when the pressure variations are lower. 

Management of the time intervals consists in connectin g connects the duration of the 
15 time interval n+1 to the maximum pressure variation observed during time interval n. 
The user must therefore provide the following data : 

Initial time interval : dt^ 
Minimum time interval : dt^ 
Maximum time interval : dt^ 
20 Lower limit of pressure variation : dp min 



10 



c) Time interval management 



Upper limit of pressure variation : dp, 
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Time interval increase factor : rdt + (>1) 
Time interval reduction factor : rdt - (<1). 

From these values, management of the time interval is performed as follows : 

At the time t„, the majorant maioritv of the pressure variations in the fracture meshes 
and the matrix blocks during time interval n is calculated (i.e. between times t 1 "' and t°). 
According to the value of this pressure variation dp, the time interval is either : 

increased by a factor rdt+ if dp < dp min , 

decreased by a factor rdt- if dp > dp max , 

unchanged in the other cases. 

After each well flow rate condition change, the time interval is re-initialized to the 
value dti„i. Furthermore, an optimization is performed at the end of the simulation 
period : if tf is the time to be simulated to reach the end of the period and dt the planned 
time interval, the time interval selected is min(tp2,dt). 

5) Simulation input data 

a) Fracture network 

The geometry of the fracture network and the attributes of the fractures 
(conductivity, opening) are given in the form of a file as in the method described in the 
aforementioned US patent FR 2.757.947 6.023.656 . 

b) Petrophysics 

The petrophysical data to be provided are as follows (the unit is given between 
parentheses) : 
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compressibility of the fracture network (c f in bar* 1 ) 
compressibility of the matrix (c m .in bar* 1 ) 
matrix permeability (K in mD) 

matrix porosity dim e n s ion l c s s ^ for computing the pore volume . 

5 The data in question are considered to be homogeneous on the fractured medium 

considered. 

c) Fluids 

The data relative to the fluids contained in the fractured medium concern the oil 
(mobile) and the water (immobile) that is always present in residual amounts in the 
10 matrix : 

residual water saturation in the matrix (Sw dimensionless) 
compressibility of the water (c w in bar 1 ) 
compressibility of the oil (c 0 in bar 1 ) 
viscosity of the oil (u^ in cpo). 

15 d) Well 

The data relative to the well are : 

its geometry, in the form of a series of connected segments 

the imposed flow rates, in the form of a curve giving the imposed flow rate as a 
function of time. 
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e) Time management 

The values to be provided are the values already defined in the chapter relative to 
the time interval management : dt^, dt,™, dt^, dp mim dp mixJ rdt+ and rdt-. 

Comparative meshing example 

The following comparative example illustrates the meshing simplification and the 
correlative saving in flow calculating time that can be obtained by implementing the 
method. By way of comparison, it has been observed that the fractured zone shown in 
Fig.4, that normally would have been meshed with 80,000 meshes with conventional 
meshing techniques, can be meshed here with about 4000 meshes only. 

Fig. 5 illustrates a flow chart of the overall processing of the invention. Processing 
starts at step 10. where the fracture network is defined as discussed above under 
fracture network definition. Processing proceeds to step 12 where the well is defined as 
discussed above under well representation. Processing proceeds to step 14 where the 
matrix medium is defined as discussed above with respect to matrix medium definition. 
Processing proceeds to step 16 where computation of the transmissivities occurs as 
discussed above under matrix-fracture transmissivitv. elsewhere with respect to the 
transmissitivity Tj and as further discussed with respect to a linear system as set forth in 
equation 9. Processing proceeds to step 18 where initialization occurs as discussed 
above with reference to fracture network and matrix medium discretization. Finally, 
processing proceeds to step 20 where a dynamic simulation loop occurs involving time 
definition. 

Fig. 6 illustrates a flow chart of the dynamic simulation loop 20. Processing in the 
dynamic simulation loop 20 begins with step 100 where building of the linear system 
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using data at time steps n. occurs as discussed above regarding time interval solution 
algorithm. Processing proceeds to step 102 where resolution of the linear system is step 
n+1 occurs as discussed under time interval solution algorithm. Processing proceeds to 
step 104 where time step management and update of well rate occurs as discussed under 
time interval management. Processing proceeds to decision point 1 06 where 
determination is made if an end of simulation has occurred. If an end of simulation is 
not sensed at decision point 106 as described above under time interval management , 
the processing returns to point 100. If the end of simulation is sensed at decision point 
106, an output of the simulation which represents the pressure evolution of the well are 
produced as indicated at step 108. 




CLAIMS 

1) A method for modelling fluid flows in a fractured multilayer porous medium in 
order to simulate interactions between pressure and flow rate variations in a well 
running across said medium, characterized in that the fractured medium is discretized 

5 by means of a mesh pattern wherein the fracture meshes are centred on nodes at the 
various fracture intersections, each node being associated with a matrix volume, and 
flows are determined between each fracture mesh and the associated matrix volume in a 
pseudosteady state. 

2) A method as claimed in claim 1, characterized in that the matrix volume 
10 associated with each fracture mesh is delimited in each layer by all of the points that are 

closer to the corresponding node than to neighbouring nodes. 

3) A method as claimed in claim 2, characterized in that each fractured layer is 
discretized in pixels and the matrix volume associated with each fracture mesh is 
delimited by determining the distance from each pixel to the closest fracture mesh. 

15 4) A method as claimed in any one of the previous claims, characterized in that the 

transmissivity value is determined for each fracture mesh - matrix block pair by 
considering that the pressure varies linearly as a function of the distance from the point 
considered to the fracture mesh associated with the block. 
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FIGURE 1 



Mai li e fissur e : fracture mesh 

Nocuds do caluul (int erse ctions) : computing nodes (intersections) 



FIGURE 2 



Maillo fi ss ur e : fracture mesh 
Bloc mutricici : matrix b i ock 

Nocuds de calcul (int e rs e ctions) : computing nodes (inters e ction s ) 



FIGURE 3 




